#include "mesh.h"
using namespace std;

Mesh::~Mesh(){}
//------------------------------------------------------------------------------
Mesh::Mesh(ClassParameter& param) {

  // --- clock info ------------------------------------------------------------
  
  double tstart = omp_get_wtime();
  
  // --- get file name ---------------------------------------------------------
  
  string fileName = param.getString("MeshFileName");
  
  cout << "[Mesh] - mesh file -------------> " << fileName << endl;	
	cout << "[Mesh] - loading mesh..." ;

  // --- Open file -------------------------------------------------------------
  
	ifstream in;
	in.open(fileName.c_str());
	if(!in) {cerr << "<<<< ERROR >>>> [Mesh] couldn't open mesh file: " << fileName;exit(1);}

  // --- Load nodes ------------------------------------------------------------
  
	string tmp;
	for(int i=0;i<10;i++){
    in >> tmp;
    if(tmp=="$Nodes"){in >> nbNodes;break;}
    if(i==9){cerr << "<<<< ERROR >>>> [Mesh] 0 nodes found in mesh file: " << fileName;exit(1);}
  }
	
  nodes = DGMatrix(nbNodes,2);
  for(int i=0;i<nbNodes;i++){in >> tmp >> nodes(i,0) >> nodes(i,1) >> tmp;} //printf("%g %g\n",nodes(i,0),nodes(i,1));}

  // --- Load elements ---------------------------------------------------------
  
  int nbTotalElmts;
  for(int i=0;i<10;i++){
    in >> tmp;
    if(tmp=="$Elements"){in >> nbTotalElmts;break;}
    if(i==9){cerr << "<<<< ERROR >>>> [Mesh] 0 elements found in mesh file: " << fileName;exit(1);}
  }
  
  DGIntMatrix* tmpETag     = new DGIntMatrix(nbTotalElmts,12);                  // temporary container
  vector<int> tmpNbElmt(3); tmpNbElmt[0]=0; tmpNbElmt[1]=0; tmpNbElmt[2]=0;     // temporary container for nbElmt: Line,Triangle,Quadrangle
  
  for(int i=0;i<nbTotalElmts;i++){
    int cpt=0;
    for(int k=0;k<3;k++,cpt++)                in >> (*tmpETag)(i,k);            // GmshTag , EType , nbTags
    for(int k=0;k<(*tmpETag)(i,2);k++,cpt++)  in >> (*tmpETag)(i,3+k);          // unused tags , geoTag
    switch ((*tmpETag)(i,1)){
      case 15:                 for(int k=cpt;k<cpt+1;k++) {in >> (*tmpETag)(i,k);}; break;
      case 1:  tmpNbElmt[0]++; for(int k=cpt;k<cpt+2;k++) {in >> (*tmpETag)(i,k);}; break;
      case 2:  tmpNbElmt[1]++; for(int k=cpt;k<cpt+3;k++) {in >> (*tmpETag)(i,k);}; break;
      case 3:  tmpNbElmt[2]++; for(int k=cpt;k<cpt+4;k++) {in >> (*tmpETag)(i,k);}; break;
      default: continue;break;
    }
  }
  
  in.close();
  
  // --- Create containers in memory -------------------------------------------
  
  int Tpos = 0;
  int Qpos = 0;
  int cpt  = 0;
  
  for(int i=1;i<3;i++){
    if(tmpNbElmt[i]>0){
      cpt++;Qpos=cpt-1;
      elmtGmshType.push_back((elmtType)(i+1));
      nbElmt.push_back(tmpNbElmt[i]);
      nbElmtFaces.push_back(i+2);
      DGIntMatrix* tmpM = new DGIntMatrix(tmpNbElmt[i],i+2);
      elmtNodesTag.push_back(*tmpM);
      vector<int> tmpEGmshTagV(tmpNbElmt[i]);
      vector<DGVector> tmpvecDGV;
      vector<DGMatrix> tmpvecDGVID;
      vector<DGMatrix> tmpvecDGMnode;
      vector<DGMatrix> tmpvecDGMJ;      
      vector<DGMatrix> tmpvecDGMinvJ; 
      vector<double>   tmpvecDetJ;
      for(int n=0;n<tmpNbElmt[i];n++){
        tmpvecDetJ.push_back(0.);
        DGMatrix tmpDGVID(3,2); tmpvecDGVID.push_back(tmpDGVID);                //NOT YET IMPLEMENTED for quad=4;
        DGVector tmpDGV(2); tmpvecDGV.push_back(tmpDGV);
        DGMatrix tmpDGMnode(i+2,2); tmpvecDGMnode.push_back(tmpDGMnode);
        DGMatrix tmpDGMinvJ(2,2); tmpvecDGMinvJ.push_back(tmpDGMinvJ);
        DGMatrix tmpDGMJ(2,2); tmpvecDGMJ.push_back(tmpDGMJ);
      } 
      elmtGmshTag.push_back(tmpEGmshTagV);
      elmtH.push_back(tmpvecDGV);
      elmtFaceID.push_back(tmpvecDGVID);
      elmtNPcoord.push_back(tmpvecDGMnode);
      elmtJ.push_back(tmpvecDGMJ);
      elmtInvJ.push_back(tmpvecDGMinvJ);
      elmtDetJ.push_back(tmpvecDetJ);
    }
  }
  
  // --- Fill containers -------------------------------------------------------

  DGIntMatrix* tmpbF = new DGIntMatrix(tmpNbElmt[0],2);
  bFacesNodesTag = *tmpbF;
  
  //cout << bFacesNodesTag.size1() << " " << bFacesNodesTag.size2() << endl;
  //cout << bFacesNodesTag(0,0) << endl;
  
  vector<int> tmpNbElmtcpt(3); tmpNbElmtcpt[0]=0; tmpNbElmtcpt[1]=0; tmpNbElmtcpt[2]=0; 
  
  for(int i=0;i<nbTotalElmts;i++){
    int EType = (*tmpETag)(i,1);
    switch (EType) {
      case 1: 
        bGmshGeoTag.push_back((*tmpETag)(i,(*tmpETag)(i,2)+2));
        for(int k=0;k<2;k++) bFacesNodesTag(tmpNbElmtcpt[0],k) = (*tmpETag)(i,(*tmpETag)(i,2)+3+k)-1;       // -1 to start node numbering from 0 and not 1 
        tmpNbElmtcpt[0]++;
        break;
      case 2: 
        elmtGmshTag[Tpos][tmpNbElmtcpt[1]] = (*tmpETag)(i,0);
        for(int k=0;k<3;k++) elmtNodesTag[Tpos](tmpNbElmtcpt[1],k) = (*tmpETag)(i,(*tmpETag)(i,2)+3+k)-1;
        tmpNbElmtcpt[1]++;
        break;
      case 3: 
        elmtGmshTag[Qpos][tmpNbElmtcpt[2]] = (*tmpETag)(i,0);
        for(int k=0;k<4;k++) elmtNodesTag[Qpos](tmpNbElmtcpt[2],k) = (*tmpETag)(i,(*tmpETag)(i,2)+3+k)-1;
        tmpNbElmtcpt[2]++;
        break;
      default: continue; break;
        break;
    }
  }
  
  // --- Delete tmp containers -------------------------------------------------
  
  delete tmpETag;
  
  // --- Set nb Total elements -------------------------------------------------
  
  nbGlobalElmts = 0;
  for(vector<int>::const_iterator eType=nbElmt.begin();eType!=nbElmt.end();++eType) nbGlobalElmts += *eType;

  // --- Set index -------------------------------------------------------------
  
  vector<int> FaceNum; FaceNum.push_back(0);
  vector<int> TypeNum; TypeNum.push_back(-1);
  
  cpt=0;
  for(vector<elmtType>::const_iterator eType=elmtGmshType.begin();eType!=elmtGmshType.end();++eType,cpt++){
    switch (*eType) {
      case Triangle:{
        DGIntMatrix* tmpI = new DGIntMatrix(3,2);
        (*tmpI)(0,0) = 0; (*tmpI)(0,1) = 1;
        (*tmpI)(1,0) = 1; (*tmpI)(1,1) = 2;
        (*tmpI)(2,0) = 2; (*tmpI)(2,1) = 0;
        nodesIndex.push_back(*tmpI);
        for(int k=0;k<3;k++){FaceNum.push_back(k);TypeNum.push_back(cpt);}
        break;}
      case Quadrangle:{
        DGIntMatrix* tmpI = new DGIntMatrix(4,2);
        (*tmpI)(0,0) = 0; (*tmpI)(0,1) = 1;
        (*tmpI)(1,0) = 1; (*tmpI)(1,1) = 2;
        (*tmpI)(2,0) = 2; (*tmpI)(2,1) = 3;  
        (*tmpI)(3,0) = 3; (*tmpI)(3,1) = 0;  
        nodesIndex.push_back(*tmpI);
        for(int k=0;k<4;k++){FaceNum.push_back(k);TypeNum.push_back(cpt);}
        break;}
      default:
        break;
    }
  }
  
  // --- Get faces connectivity ------------------------------------------------
  
  nbIntFaces   = 0;
  nbBFaces     = 0;
  DGIntMatrix* tmpbFtag = new DGIntMatrix(tmpNbElmt[0],4);
  bFacesTag    = *tmpbFtag;
  for(int i=0;i<tmpNbElmt[0];i++) bFacesNormal.push_back(DGVector(3));

  cpt = 0;
  vector<ANNpointArray> dataFaceTag;
  vector<vector<bool> > checkElmt;
  vector<ANNkd_tree*>   kdTree;

  int nbdataFace=0;
  if(tmpNbElmt[0]>0){
    dataFaceTag.push_back(annAllocPts(tmpNbElmt[0],2));
    vector<bool> tmpcheck;
    for(int i=0;i<tmpNbElmt[0];i++){
      ANNcoord* atmp = new ANNcoord[2];
      atmp[0] = bFacesNodesTag(i,0);
      atmp[1] = bFacesNodesTag(i,1);
      dataFaceTag[nbdataFace][i] = atmp;
      tmpcheck.push_back(false);
    }
    ANNkd_tree* tmpkdTree = new ANNkd_tree(dataFaceTag[nbdataFace],tmpNbElmt[0],2);
    kdTree.push_back(tmpkdTree);
    checkElmt.push_back(tmpcheck);
    nbdataFace++;
  }
  
  int eT=0;
  for(vector<elmtType>::const_iterator eTIt=elmtGmshType.begin();eTIt!=elmtGmshType.end();++eTIt,eT++){
    for(int n=0;n<nbElmtFaces[eT];n++){//nbface
      dataFaceTag.push_back(annAllocPts(nbElmt[eT],2));
      vector<bool> tmpcheck;
      for(int i=0;i<nbElmt[eT];i++){
        ANNcoord* atmp = new ANNcoord[2];
        atmp[0] = elmtNodesTag[eT](i,nodesIndex[eT](n,0));    
        atmp[1] = elmtNodesTag[eT](i,nodesIndex[eT](n,1));
        dataFaceTag[nbdataFace][i] = atmp;
        tmpcheck.push_back(false);
      }
      ANNkd_tree* tmpkdTree = new ANNkd_tree(dataFaceTag[nbdataFace],nbElmt[eT],2);
      kdTree.push_back(tmpkdTree);
      checkElmt.push_back(tmpcheck);
      nbdataFace++;
    }
  }
  
  ANNidxArray  nnIdx       = new ANNidx[1];
  ANNdistArray dists       = new ANNdist[1];
  ANNpoint cibleIntNodeTag = new ANNcoord[2];
  ANNpoint cibleBCNodeTag  = new ANNcoord[2];
  ANNdist  d               = 1e-6;
  
  DGIntMatrix* tmpITag    = new DGIntMatrix(nbGlobalElmts*4,6);
  DGMatrix*    tmpINormal = new DGMatrix(nbGlobalElmts*4,3);
  
  vector<vector<bool> >::const_iterator elmtTypeIter = checkElmt.begin();
  
  int nT=0;
  for(;elmtTypeIter!=checkElmt.end();++elmtTypeIter,nT++) {
    
    if(nT==0) continue; //skip Line
    
    vector<bool>::const_iterator elmtIter = elmtTypeIter->begin();
    
    int nE=0;
    for(;elmtIter!=elmtTypeIter->end();++elmtIter,nE++) {
      
      if(checkElmt[nT][nE]==true) continue;
      
      cibleBCNodeTag[0]  = dataFaceTag[nT][nE][0];
      cibleBCNodeTag[1]  = dataFaceTag[nT][nE][1];
      cibleIntNodeTag[0] = dataFaceTag[nT][nE][1];
      cibleIntNodeTag[1] = dataFaceTag[nT][nE][0];
      
      for(int k=1;k<nbdataFace;k++) {
        nnIdx[0]= -1;
        kdTree[k]->annkFRSearch(cibleIntNodeTag,d,1,nnIdx,dists,0);

        if(nnIdx[0]!=-1){
          (*tmpITag)(nbIntFaces,0) = TypeNum[nT]; 
          (*tmpITag)(nbIntFaces,1) = nE;
          (*tmpITag)(nbIntFaces,2) = FaceNum[nT];
          (*tmpITag)(nbIntFaces,3) = TypeNum[k];
          (*tmpITag)(nbIntFaces,4) = nnIdx[0];
          (*tmpITag)(nbIntFaces,5) = FaceNum[k];
          (*tmpINormal)(nbIntFaces,0) = -nodes(elmtNodesTag[TypeNum[nT]](nE,nodesIndex[TypeNum[nT]](FaceNum[nT],0)),1)+nodes(elmtNodesTag[TypeNum[nT]](nE,nodesIndex[TypeNum[nT]](FaceNum[nT],1)),1);
          (*tmpINormal)(nbIntFaces,1) = nodes(elmtNodesTag[TypeNum[nT]](nE,nodesIndex[TypeNum[nT]](FaceNum[nT],0)),0)-nodes(elmtNodesTag[TypeNum[nT]](nE,nodesIndex[TypeNum[nT]](FaceNum[nT],1)),0);
          (*tmpINormal)(nbIntFaces,2) = sqrt((*tmpINormal)(nbIntFaces,0)*(*tmpINormal)(nbIntFaces,0)+(*tmpINormal)(nbIntFaces,1)*(*tmpINormal)(nbIntFaces,1));
          (*tmpINormal)(nbIntFaces,0) = (*tmpINormal)(nbIntFaces,0)/(*tmpINormal)(nbIntFaces,2);
          (*tmpINormal)(nbIntFaces,1) = (*tmpINormal)(nbIntFaces,1)/(*tmpINormal)(nbIntFaces,2);
          checkElmt[nT][nE]=true;
          checkElmt[k][nnIdx[0]]=true;
          nnIdx[0]=-1;
          nbIntFaces++;
          break;
        }else if(k==nbdataFace-1){
          nnIdx[0]= -1;
          kdTree[0]->annkFRSearch(cibleBCNodeTag,d,1,nnIdx,dists,0);
          if(nnIdx[0]!=-1){
            bFacesTag(nbBFaces,0) = TypeNum[nT]; 
            bFacesTag(nbBFaces,1) = nE;
            bFacesTag(nbBFaces,2) = FaceNum[nT];
            bFacesTag(nbBFaces,3) = bGmshGeoTag[nnIdx[0]];
            bFacesNormal[nbBFaces](0) = -nodes(elmtNodesTag[TypeNum[nT]](nE,nodesIndex[TypeNum[nT]](FaceNum[nT],0)),1)+nodes(elmtNodesTag[TypeNum[nT]](nE,nodesIndex[TypeNum[nT]](FaceNum[nT],1)),1);
            bFacesNormal[nbBFaces](1) = nodes(elmtNodesTag[TypeNum[nT]](nE,nodesIndex[TypeNum[nT]](FaceNum[nT],0)),0)-nodes(elmtNodesTag[TypeNum[nT]](nE,nodesIndex[TypeNum[nT]](FaceNum[nT],1)),0);
            bFacesNormal[nbBFaces](2) = sqrt(bFacesNormal[nbBFaces](0)*bFacesNormal[nbBFaces](0)+bFacesNormal[nbBFaces](1)*bFacesNormal[nbBFaces](1));
            bFacesNormal[nbBFaces](0) = bFacesNormal[nbBFaces](0)/bFacesNormal[nbBFaces](2);
            bFacesNormal[nbBFaces](1) = bFacesNormal[nbBFaces](1)/bFacesNormal[nbBFaces](2);
            checkElmt[nT][nE]=true;
            checkElmt[0][nnIdx[0]]=true;
            nnIdx[0]=-1;
            nbBFaces++;
            break;
          }
        }
      
      }
    }
    
  }
  
  DGIntMatrix* tmpIFTag = new DGIntMatrix(nbIntFaces,6);
  intFacesTag    = *tmpIFTag;
  for(int i=0;i<nbIntFaces;i++) intFacesNormal.push_back(DGVector(3));

  for(int i=0;i<nbIntFaces;i++){
    for(int j=0;j<6;j++){
      intFacesTag(i,j) = (*tmpITag)(i,j);
      if(j<3) intFacesNormal[i](j) = (*tmpINormal)(i,j);
    }
  }
  //for bFacesNormal it's already computed in kdree loop
  
  // --- add fIDs for each element ---------------------------------------------
  
  for(int fID=0;fID<nbIntFaces;fID++){
    
    int eTL      = getIntFacesTag(fID,0);                               
    int eIDL     = getIntFacesTag(fID,1); 
    int eFNumL   = getIntFacesTag(fID,2);
    int eTR      = getIntFacesTag(fID,3);                               
    int eIDR     = getIntFacesTag(fID,4); 
    int eFNumR   = getIntFacesTag(fID,5);

    elmtFaceID[eTL][eIDL](eFNumL,0) = fID;
    elmtFaceID[eTL][eIDL](eFNumL,1) = 0;
    elmtFaceID[eTR][eIDR](eFNumR,0) = fID;
    elmtFaceID[eTR][eIDR](eFNumR,1) = 1;
    
  }

  for(int bID=0;bID<nbBFaces;bID++){
    
    int eTL      = getBFacesTag(bID,0);                               
    int eIDL     = getBFacesTag(bID,1); 
    int eFNumL   = getBFacesTag(bID,2);
    
    elmtFaceID[eTL][eIDL](eFNumL,0) = bID;
    elmtFaceID[eTL][eIDL](eFNumL,1) = 3;
    
  }
  
  // --- Check if connectivity is correctly done ------------------------------- // NOT YET IMPLEMENTED for quadrangle

  bool passline = false;
  bool checkVal = false;
  nT=0;
  for(vector<vector<bool> >::const_iterator fNumIt=checkElmt.begin()+1;fNumIt!=checkElmt.end();++fNumIt,nT++){
    if(passline==false) {passline=true;continue;} 
    int nE=0;
    for(vector<bool>::const_iterator eId=fNumIt->begin();eId!=fNumIt->end();++eId,nE++){
      if(*eId==false){
        cout << "Face = " << FaceNum[nT] << "   GmshTag = " << elmtGmshTag[0][nE] << " nE--->" << nE << endl;
        checkVal=true;
      }
    }
  }
  if(checkVal==true){cerr << endl << "<<<< ERROR >>>> [Mesh] some faces have no connectivity" << endl;exit(1);}

  /*
  DGVector checkNbF(nbElmt[0]);                                                 
  
  for(int i=0;i<nbIntFaces;i++){
    checkNbF(intFacesTag(i,1)) += 1; 
    checkNbF(intFacesTag(i,4)) += 1;
  }
  for(int i=0;i<nbBFaces;i++){
    checkNbF(bFacesTag(i,1)) += 1;
  }
  for(int i=0;i<nbElmt[0];i++){
    if(checkNbF(i)!=3)
    cout << checkNbF(i) << "-->" << elmtGmshTag[0][i] << endl;
    //if(checkNbF(i)!=3) {cerr << endl << "<<<< ERROR >>>> [Mesh] face connectivity is <3 for at least one element!" << endl;exit(1);}
  }
  */

  // --- Add periodicity -------------------------------------------------------
  // not yet implemented

  // --- memory free -----------------------------------------------------------
  
  delete tmpITag;
  delete tmpINormal;
  for(int i=0;i<nbdataFace;i++) delete(kdTree[i]);
  
  // --- Build Jacobian matrices -----------------------------------------------
  
  buildJacobian();
  
  // --- Build hmax & hmin -----------------------------------------------------
  
  computeElmtSize();
  
  // --- infos -----------------------------------------------------------------
  
  double tend = omp_get_wtime();
  double elapsedTime = tend-tstart;
  
  int order = param.getInteger("Order");
  int nbEffDof = 0;
  switch (order) {
    case 1: nbEffDof = nbNodes;break;
    case 2: nbEffDof = nbNodes + nbBFaces + nbIntFaces;break;
    case 3: nbEffDof = nbNodes + 2*(nbBFaces+nbIntFaces) + nbGlobalElmts;break;      
  }
	cout << "\r[Mesh] - nb of nodes -----------> " << nbNodes << endl;
  cout << "[Mesh] - nb of elements --------> " << nbGlobalElmts << "   (" << tmpNbElmt[1] << " triangles," << tmpNbElmt[2] << " quadrangles)" << endl;
  cout << "[Mesh] - nb of interfaces ------> " << nbIntFaces << endl;
  cout << "[Mesh] - nb of BC faces --------> " << nbBFaces << endl;
  cout << "[Mesh] - nb DOF ----------------> " << 0.5*(order+1)*(order+2)*nbGlobalElmts << endl;
  cout << "[Mesh] - nb eff. DOF -----------> " << nbEffDof << "   (you have to subsract the number of geom nodes!) " << endl;
  cout << "[Mesh] - inscribed circle ------> (" << Hmaxmin(0,0) << "," << Hmaxmin(0,1) << ")     rmin,rmax" << endl;
  cout << "[Mesh] - circumscribed circle --> (" << Hmaxmin(1,0) << "," << Hmaxmin(1,1) << ")     rmin,rmax" << endl; 
  cout << "[Mesh] - loading mesh ----------> done in " << elapsedTime << " seconds"<< endl ;
  
}
//------------------------------------------------------------------------------
void Mesh::buildJacobian(){
  
  int eT=0;
  for(vector<vector<DGMatrix> >::iterator eTIt=elmtJ.begin();eTIt!=elmtJ.end();++eTIt,eT++){
    
    int nbNP    = elmtNPcoord[eT][0].size1(); 
    
    int eID=0;
    for(vector<DGMatrix>::iterator eJIt=eTIt->begin();eJIt!=eTIt->end();++eJIt,eID++){
      
      for(int i=0;i<nbNP;i++) DGVector(elmtNPcoord[eT][eID].touchRow(i)).memcpy(DGVector(nodes.seeRow(elmtNodesTag[eT](eID,i))));
      
      for(int i=0;i<2;i++){
        for(int j=0;j<2;j++){
          for(int k=0;k<nbNP;k++){
            elmtJ[eT][eID](i,j) += nodes(elmtNodesTag[eT](eID,k),j)*dus(eT,k,i);
          }
        }
      }
      
      elmtInvJ[eT][eID].memcpy(elmtJ[eT][eID]);
      elmtInvJ[eT][eID].invert();
      elmtDetJ[eT][eID] = 0.5*(elmtJ[eT][eID](0,0)*elmtJ[eT][eID](1,1)
                         - elmtJ[eT][eID](0,1)*elmtJ[eT][eID](1,0));              // detJ = Area for a triangle
    }
  } 
      
}
//------------------------------------------------------------------------------
double Mesh::dus(int eType,int r,int q){
  
  switch (elmtGmshType[eType]) {
    case Triangle:{
      double coef[3][3] = {{1,-1,-1},{0,1,0},{0,0,1}};
      double tmp[3][2]  = {{0,0},{0,0},{0,0}};
      for(int i=0;i<3;i++){
        for(int j=0;j<3;j++){
          tmp[i][0] += coef[i][j]*dmonomeof(eType,0.1,0.1,j,0); 
          tmp[i][1] += coef[i][j]*dmonomeof(eType,0.1,0.1,j,1);
        }
      }
      return tmp[r][q];
      break;}
    case Quadrangle:{
      //NOT YET IMPLEMENTED
      break;}
  }
  
  return 0.;
  
}
//------------------------------------------------------------------------------
double Mesh::monomeof(int eType,double u,double v,int i){
  
  switch (elmtGmshType[eType]) {
    case Triangle:{
      double pp[3][2] = {{0,0},{1,0},{0,1}};
      return pow(u,pp[i][0])*pow(v,pp[i][1]);}
    case Quadrangle:{
      //NOT YET IMPLEMENTED
      break;}
  }
  
  return 0.;
  
}
//------------------------------------------------------------------------------
double Mesh::dmonomeof(int eType,double u,double v,int i,int k){
  
  switch (elmtGmshType[eType]) {
    case Triangle:{
      double pp[3][2] = {{0,0},{1,0},{0,1}};
      if(pp[i][k]<=0){
        return 0.;
      }else{
        switch (k) {
          case 0:
            return pp[i][k]*pow(u,pp[i][0]-1.)*pow(v,pp[i][1]);
          case 1:
            return pp[i][k]*pow(u,pp[i][0])*pow(v,pp[i][1]-1.);
        }
      }
      break;}
    case Quadrangle:{
      //NOT YET IMPLEMENTED
      break;}
  }
  
  return 0.;

}
//------------------------------------------------------------------------------
void Mesh::computeElmtSize(){

  //NOT YET IMPLEMENTED for quadrangle
  
  Hmaxmin = DGMatrix(2,2);
  Hmaxmin(0,0) = 1000.; //hmin-min
  Hmaxmin(0,1) = 0.;    //hmin-max
  Hmaxmin(1,0) = 1000.; //hmax-min
  Hmaxmin(1,1) = 0.;    //hmax-max
  
  int eT=0;
  for(vector<vector<DGVector> >::iterator eTIt=elmtH.begin();eTIt!=elmtH.end();++eTIt,eT++){
    
    int nbNP    = elmtNPcoord[eT][0].size1(); 
    
    int eID=0;
    for(vector<DGVector>::iterator eHIt=eTIt->begin();eHIt!=eTIt->end();++eHIt,eID++){
      
      double    detJ    = getDetJ(eT,eID);
      DGMatrix* NPcoord = getNPcoord(eT,eID);
      DGVector  GPcoord(2);
      for(int i=0;i<nbNP;i++) for(int dim=0;dim<2;dim++) GPcoord(dim) += elmtNPcoord[eT][eID](i,dim);
      for(int dim=0;dim<2;dim++) GPcoord(dim) /= nbNP;
      
      double perimeter = 0.;
      for(int i=0;i<nbNP;i++){
        perimeter += sqrt(pow(elmtNPcoord[eT][eID](nodesIndex[eT](i,1),0)-elmtNPcoord[eT][eID](nodesIndex[eT](i,0),0),2)
                         +pow(elmtNPcoord[eT][eID](nodesIndex[eT](i,1),1)-elmtNPcoord[eT][eID](nodesIndex[eT](i,0),1),2));
      }
      
      elmtH[eT][eID](0) = 2*detJ/perimeter;                                     // hmin --> inscribed circle r=2*A/P    or A=detJ
      elmtH[eT][eID](1) = sqrt(pow(GPcoord(0)-elmtNPcoord[eT][eID](0,0),2)
                              +pow(GPcoord(1)-elmtNPcoord[eT][eID](0,1),2));    // hmax --> circumscribed circle
      
      if(elmtH[eT][eID](0)<Hmaxmin(0,0)) Hmaxmin(0,0)=elmtH[eT][eID](0);
      if(elmtH[eT][eID](0)>Hmaxmin(0,1)) Hmaxmin(0,1)=elmtH[eT][eID](0);
      if(elmtH[eT][eID](1)<Hmaxmin(1,0)) Hmaxmin(1,0)=elmtH[eT][eID](1);
      if(elmtH[eT][eID](1)>Hmaxmin(1,1)) Hmaxmin(1,1)=elmtH[eT][eID](1);
            
    }
  }

}